# ============================================================
# Thomas König & Stefan Eschenwecker
# The European Court of Justice and legal European integration
# 
# This script produces the regression model object 
# for the Combined Model (see Table 1 in main text)
# and Figure A4 in online appendix.
# ============================================================

# load packages
library(rstan)
library(brms)
library(tidyverse)

# set options for the script
options(mc.cores = parallel::detectCores())
options(scipen = 10)
theme_set(theme_bw(base_size = 14))


# load data
CourtData <- read.csv("Court/Data/PreparedCourtData.csv")

# prepare data
# (DV + IVs as ordered factors)
# (remove cases where all actors claimed not applicable)
CourtData <- CourtData %>%
  mutate(
    pos_cjeu_ln = factor(pos_cjeu_ln,
      ordered = T,
      levels = c("PS", "Ambi", "ME")
    ),
    supra_signal = factor(supra_signal,
      ordered = F,
      levels = c( "Uninfo", "Cont", "Weak PS", "Strong PS", "Weak ME", "Strong ME")
    ),
    supra_signal_int = factor(supra_signal_int,
      ordered = T,
      levels = c("No Intensity", "Weak Intensity", "Strong Intensity"))
  ) %>%
  filter(all_not_applic == 0)


# sample priors
# (Hint: takes some time to run)
PriorModel <- brm(
  bf( # no intercept for model identification
    pos_cjeu_ln ~ supra_signal + net_ms_pref_std +
      share_judges_std + (1 | iuropa_case_id)
  ) +
    lf(
      disc ~ 0 + mo(supra_signal_int) + ms_signal_int_std +
      share_judges_std + (1 | iuropa_case_id),
       cmc = F), # no intercept for model identification
  family = cumulative(link = "probit"),
  prior = c(
    # default alpha = 1 Dirichlet prior for monotonic effect
    prior(normal(0, 1), class = b), # standardized location predictors
    prior(exponential(1 / 0.463), class = sd, group = iuropa_case_id), # 95% interval of 0.25 to 4
    prior(normal(-0.43, 2), class = Intercept, coef = 1), # location is qnorm of 1/3
    prior(normal(0.43, 2), class = Intercept, coef = 2), # location is qnorm of 2/3
    prior(normal(0, log(4) / 2), class = b, dpar = disc), # expect difference by factor of max 4
    prior(exponential(1 / 0.463), class = sd, dpar = disc) # 95% interval of 0.25 to 4
    ),
  init = 0, # avoid rejected starting values
  data = CourtData,
  sample_prior = "only",
  chains = 4, cores = 4, warmup = 1000, iter = 2500,
  seed = 1408, control = list(adapt_delta = 0.99),
  backend = "cmdstanr"
)


# plot prior predictive distribution
# (Figure A4 in online appendix)
set.seed(1408)
PriorPredPlot <- pp_check(PriorModel,
  type = "bars",
  ndraws = 1000, prob = 0.95
) +
  ylab("Count") +
  xlab("Outcome") +
  scale_x_continuous(
    labels = c("Preserve Sovereignty", "Ambivalent", "More Europe"),
    breaks = c(1, 2, 3)
  ) +
  scale_colour_manual(values = "grey20") +
  scale_fill_manual(values = "grey80") +
  theme(legend.position = "none")
PriorPredPlot$layers[[1]]$aes_params$colour <- "grey40"
PriorPredPlot$layers[[2]]$aes_params$linewidth <- 1.2
PriorPredPlot$layers[[2]]$aes_params$size <- 1.2
PriorPredPlot

# save plot
ggsave(plot = PriorPredPlot, device = "pdf",
       filename = "PriorPredictiveDistribution.pdf",
       path = "Court/Plots/", height = 6, width = 8)



# estimate regression model
# (automatically save the object)
# (Hint: takes some time to run)
OrdProbModel <- brm(
  bf(
    pos_cjeu_ln ~ supra_signal + net_ms_pref_std +
      share_judges_std + (1 | iuropa_case_id)
  ) +
    lf(
      disc ~ 0 + mo(supra_signal_int) + ms_signal_int_std +
      share_judges_std + (1 | iuropa_case_id),
       cmc = F),
  family = cumulative(link = "probit"),
  prior = c(
    # default alpha = 1 Dirichlet prior for monotonic effect
    prior(normal(0, 1), class = b), # standardized location predictors
    prior(exponential(1 / 0.463), class = sd, group = iuropa_case_id), # 95% interval of 0.25 to 4
    prior(normal(-0.43, 2), class = Intercept, coef = 1), # location is qnorm of 1/3
    prior(normal(0.43, 2), class = Intercept, coef = 2), # location is qnorm of 2/3
    prior(normal(0, log(4) / 2), class = b, dpar = disc), # expect difference by factor of max 4
    prior(exponential(1 / 0.463), class = sd, dpar = disc) # 95% interval of 0.25 to 4
  ),
  init = 0, # avoid rejected starting values
  data = CourtData,
  chains = 4, cores = 4, warmup = 1000, iter = 2500,
  seed = 1408, control = list(adapt_delta = 0.99),
  file = "Court/Models/UnequalPrecision",
  backend = "cmdstanr"
)

# print summary
summary(OrdProbModel)

# check for non-convergence
shinystan::launch_shinystan(OrdProbModel)
